Bilayer honeycomb lattices (BHLs) with Bernal stacking have been constantly studied in condensed matter physics in relation with the fabrication of single- and multi-layer graphene.^[1–4] Very recently, unconventional superconductivity with T_C ∼ 1.7 K has been observed in twisted bilayer graphene superlattices close to the first magic angle,^[5] which reboots enormous research interests upon this class of systems. Multiple two-dimensional (2D) materials with great current research interest can be described by the BHLs, e.g., the hexagonal boron nitride bilayer (h-BN),^[6] the commensurate h-BN graphene bilayer,^[7] and more complicated materials such as the Mn⁴⁺ ions in Bi₃ Mn₄O₁₂(NO₃).^[8] In a Bernal-stacked BHL, the density of states (DOS) and the local density of states (LDOS) are largely affected by the inter-layer hopping energies.^[9,10] The LDOS on the two sublattices in the same plane becomes distinct, and this inequivalence is expected to induce interesting physics in such systems.

The Kondo effect takes place when a magnetic impurity forms a singlet with the conduction electrons at a temperature lower than the Kondo temperature, and has been constantly studied by using various methods.^[11–19] The Kondo effect and the RKKY interactions in two- or three-dimensional (3D) Dirac systems have been studied intensively in recent years.^[20–24] At half-filling due to the vanishing DOS, the problem of magnetic impurity in Dirac systems falls into the category of the pseudo-gap Kondo problem.^[25–27] There exists a critical value of hybridization for the impurity and the conduction electrons to form a bound state.^[28,29]

The Kondo problem in the bilayer graphene (BLG) has been studied using mean-field theory,^[30–32] and it has been reported that the local moment of the impurity is tunable in a wide range, by varying parameters such as the inter-layer hopping energies, the chemical potential, and the external electric field. The effect of spatial inhomogeneity on the local moment formation in a bilayer graphene with only the inter-layer hopping energy connecting pairs of the dimer sites has been studied.^[31] The non-dimer inter-layer hopping energy has slightly smaller magnitude than the dimer ones but is often ignored in the theoretical studies. These non-dimer couplings can induce trigonal wrapping in the bilayer system.^[9] In comparison with BLG, BHL is a more general theoretical model in the sense that the system spans a wider range of inter-layer hopping energies covering a large class of materials with great current research interest. However, systematic studies on the local moment formation in BHLs using non-perturbative numerical methods are still lacking.

The purpose of this paper is to investigate the properties of a spin-1/2 magnetic impurity in a Bernal-stacked BHL. The host system is described by a tight-binding model with the two most significant inter-layer hoppings, t₁ between pairs of dimer sites and t₃ between pairs of non-dimer sites, which is used to describe the kinetic energy of the system. The local moment of a magnetic impurity in such systems is expected to have an asymmetric response, e.g., the local moment shows an asymmetric response to the inter-layer hoppings depending on which sublattice the impurity is coupled to.

The Hirsch–Fye quantum Monte Carlo (QMC) numerical method^[33,34] we utilize can offer non-perturbative results for interacting impurity systems. The results we present are unbiased and numerically exact, namely, they are free of systematic errors, such that the local moment can be considered quantitatively in a wide parameter range. The Hirsch–Fye QMC method has already been used successfully to study the magnetic impurity problem in normal metals,^[33–35] single layer graphene,^[36] dilute magnetic semiconductors,^[37] correlated fermions,^[38] and systems with spin–orbit coupling such as topological insulators.^[39]

We use the Anderson impurity model to investigate the local moment formation of a spin-1/2 magnetic impurity in a Bernal-stacked BHL. The model Hamiltonian consists of three parts: the kinetic energy term H₀ describing the Bernal-stacked BHL, the magnetic impurity part H_d, and the hybridization interaction H_V. The total Hamiltonian is given by 1 The tight-binding lattice model Hamiltonian of a BHL with Bernal stacking reads 2 where μ is the chemical potential of the system, and μ = 0 corresponds to the half-filled case. a_lσ(i) (b_lσ(j)) annihilates an electron with spin-σ on the i-th (j-th) site in sublattice A (B), l = 1,2 is the layer index, and ⟨ … ⟩ denotes an in-plane or inter-layer nearest neighbor (NN) bond. t is the in-plane NN hopping, and is used as the energy unit in our calculation. If both t₁ and t₃ vanish, the Hamiltonian describes two isolated graphene sheets. t₁ is the inter-layer NN hopping energy between pairs of dimer sites in sublattices A₁ and A₂ as shown in Fig. 1, which will induce the parabolic dispersion relation and finite DOS at half-filling.^[10]t₃ is the inter-layer hopping energy between pairs of non-dimer sites in sublattices B₁ and B₂, and this term induces trigonal wrapping in the bilayer system.^[9]

Fig. 1.

Figure Option ViewDownloadNew Window

Fig. 1. (color online) Schematics of a BHL with Bernal stacking. (a) t is the in-plane nearest-neighbor hopping energy, and is used as the energy unit. A1 (A2) is the sublattice containing the dimer sites with hopping energy t1 in layer 1 (layer 2), while B1 (B2) is the sublattice containing the non-dimer sites with trigonal hopping energy t3 in layer 1 (layer 2). “Case A” represents that the magnetic impurity is coupled to a site on sublattice A1, while “case B” represents the coupling with a site on sublattice B1. (b) Top view of the Bernal-stacked BHL.

The impurity part Hamiltonian is given by 3 where d_σ annihilates an electron with spin-σ at the localized states, and U is the on-site Coulomb repulsion. The impurity energy level ε_d is chosen as −U/2, which will give well-developed local moments at low temperatures.

The hybridization interaction between the magnetic impurity and the conduction electrons is given by 4 where V is the strength of the hybridization interaction. We assume that the magnetic impurity is coupled to layer 1, and the location of the impurity is denoted as i = 0. As given in Fig. 1, A₁ (A₂) is the sublattice containing the dimer sites with hopping energy t₁ in layer 1 (2), while B₁ (B₂) is the sublattice containing the non-dimer sites with trigonal hopping energy t₃ in layer 1 (2). Hereinafter, we use the term “case A” to indicate the case that the magnetic impurity is coupled to a site on sublattice A₁, and “case B” denotes the coupling with a site on sublattice B₁. Hence for “Case A” c_1σ(i = 0) = a_{1 σ}(i = 0), while for “Case B” c_1σ(i = 0) = b_1σ(i = 0) in Eq. (4).

The Hirsch–Fye algorithm we utilize naturally returns the imaginary-time Green’s function G_d(τ) = ∑_σG_dσ (τ) of the impurity. With Green function, we can easily calculate the expectation value of the local moment squared 5 based on the fact that the impurity charge can either be zero or one for a given spin. ⟨n_d⟩ = ⟨n_d↑ + n_d↓⟩ is the total charge on the impurity site, and ⟨n_d↑n_d↓⟩ describes the expectation value of the double occupancy. The closer this value is to one, the more fully developed the moment is.

Using the imaginary-time Green’s function obtained from the QMC method, we can calculate the spectral density A(ω) = ∑_σA_σ(ω) by numerically solving^[40] 6

We also use an extension of the Hirsch–Fye algorithm^[41] to compute the charge–charge correlation function 7 and the spin–spin correlation function 8 where m_i is the magnetization on the BHL with index i.

In this section, we investigate the asymmetric response of a spin-1/2 magnetic impurity to the inter-layer hopping energies in a Bernal-stacked BHL. There are several tunable parameters in the model Hamiltonian, such as the hopping energies t₁ and t₃, the on-site Coulomb repulsion U, and the hybridization strength V. In single layer graphene, the Coulomb repulsion U and the hopping energy are shown to be tunable by experimental means such as strain.^[42] In the theoretical studies on BLG, t₁ and t₃ are usually considered as the most significant inter-layer hopping terms, and the values are taken as t₁ ∼ t₃ ∼ 0.1t.^[9] Here in this present paper, we vary the inter-layer hoppings in a wider regime: both t₁ and t₃ take the values from 0 to 0.8t, in order to simulate the case that the magnetic impurity is coupled to a BHL beyond BLG. We also tune the values of U and V, in order to systematically investigate the effects of the on-site Coulomb repulsion and the hybridization on the local moment formation.

Noting that the Hirsch–Fye QMC method is based on a finite-temperature algorithm, we choose the value β = T⁻¹ = 40t⁻¹ in all the results that we show below. In order to ensure the validity of our results, we also tested with lower temperatures, and the results show minor differences. We believe that the temperatures we choose are low enough to offer correct results for the moderate-U regime we use in our simulations. The impurity energy level is always chosen as ε_d = −U/2, which will give rise to a well-developed local moment. It is well known that according to the mean-field solution of a magnetic impurity in a normal metal, the magnitude of the local moment is affected by the relative position of ε_d with respect to the Fermi level, and is also influenced by the value of U. We use “case A” to represent the case that the magnetic impurity is coupled to the site on sublattice A₁ (dimer sites), and use “case B” to denote the case that the magnetic impurity is coupled to the site on sublattice B₁ (non-dimer sites). Due to the inter-layer hopping terms t₁ and t₃, the local moment formed on the magnetic impurity for the two cases are expected to show distinctions.

In Fig. 2, we show the QMC results of versus μ for various values of t₁ and t₃. The parameters are chosen as U = 1.6t, β = 40t, V = 1.0t, and ε_d = −U/2. If both t₁ and t₃ vanish, the host material given in Eq. (2) corresponds to two decoupled single layer honeycomb lattices with linear dispersion, and μ = 0 is the charge neutral point. As μ is lowered from μ = 0, for all the cases we can see that the squared local moment decreases as well. Note that the impurity energy level ε_d is determined by the choice of U, and does not change together with μ. Hence as the chemical potential is lowered, both the impurity charge and the double occupancy are lowered. However, decreases faster, such that the local moment is tuned to a smaller value with the decreasing μ according to Eq. (5). The asymmetric response of the magnetic impurity local moment to the inter-layer hopping energies can be observed in Fig. 2. The results for “case A” are given in Fig. 2(a). The most interesting thing is the effects of t₁ and t₃ when |μ| ∼ 0. The changes of magnetic moment with respect to the inter-layer hoppings t₁ or t₃ for μ = 0 in “case A” are shown in the insets. The results for other μ values are qualitatively the same, so they are not shown. We can see that grows monotonically as t₁ increases, but decreases as t₃ increases. In a realistic system such as BLG, both t₁ and t₃ take finite values, and the magnitude of the local moment is determined by the competition of t₁ and t₃. We see that the results for t₁ = 0 and t₃ = 0.2 are almost identical to the results for t₁ = 0 and t₃ = 0. This is because the inter-layer hopping energy t₃ only connects the non-dimer sites, so for small t₃ the effects on the local moment for “case A” are not significant, although large t₃ can still decrease the local moment. In Fig. 2(b), we present the results of for “case B”. We can see that the tendency of with decreasing μ is the same as Fig. 2(a), but the magnitude of local moment for the same t₁ and t₃ is distinct from those given in Fig. 2(a), showing an asymmetric response. The dependence of on t₁ as well as t₃ is completely opposite to those given in Fig. 2(a). The squared local moment around μ = 0 decreases as t₁ increases, but grows monotonically as t₃ increases. The inset in Fig. 3(b) shows the changes of magnetic moment with respect to the inter-layer hoppings t₁ or t₃ for μ = 0 in “case B”. We can see that grows monotonically as t₃ increases, but decreases as t₁ increases. The effect of t₁ and t₃ is completely opposite to that shown in the inset of Fig. 3(a). This phenomenon is related to the changes of LDOS with respect to the inter-layer hopping energies, and we will address this issue in Section 4 with the results of spectral densities.

Fig. 2.

Figure Option ViewDownloadNew Window

Fig. 2. (color online) The results of versus μ for various values of t1 and t3. The parameters are chosen as U = 1.6t, β = 40t−1, V = 1.0t, and εd = −U/2. (a) For the “case A”, grows monotonically as t1 increases, but decreases as t3 increases in a wide range of μ. Inset shows the values of with respect to t1 or t3 for μ = 0. (b) For the “case B”, decreases as t1 increases, and grows monotonically as t3 increases. Inset shows the values of with respect to t1 or t3 for μ = 0. The inter-layer hoppings have the opposite effect on the local moment in comparison with “case A”. For all the cases, decreases as μ is lowered.

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. (color online) versus μ for various values of U. The parameters are chosen as β = 40t, V = 1.0t, t1 = 0.2t, t3 = 0.1t, and εd −μ = −U/2. The lines from top to bottom are listed in the order of decreasing U. Circles represent the results for “case A” and triangles for “case B”. Larger local moment is favored for larger U, but the relative size of for “case A” and “case B” is determined by the choice of t1 and t3.

In general, the results we give in Fig. 2 for vanishing t₃ are in consistence with the mean-field results given by previous studies,^[31] in which only the t₁ type inter-layer hopping energy is considered. The inter-layer hoppings induce inhomogeneity on the dimer and non-dimer sites, and the local moment is more easily formed on the dimer sites.

In order to study the effect of on-site Coulomb repulsion U on the local moment of the magnetic impurity, we show the results of versus μ for various values of U in Fig. 3. The parameters are chosen as β = 40t⁻¹, V = 1.0t, t₁ = 0.2t, t₃ = 0.1t, and ε_d = −U/2. The choice of t₁ and t₃ values are close to the values in the real BLG system. If we choose different values of inter-layer hopping energies, the magnitude of the local moment shall change accordingly, but the qualitative results will remain the same.

The lines from top to bottom are listed in the order of decreasing U. Circles represent the results for “case A”, and triangles for “case B”. For all the U values of concern, the magnitude of the decreases as the chemical potential is switched away from the zero point. We can see that as the Coulomb repulsion U decreases, for both the “case A” and “case B” decrease as well. This is reasonable since larger U can better exclude the double occupancy, so it favors a larger local moment. Around μ∼ 0, “case A” always presents a larger local moment in comparison with “case B” for the same set of parameters. Note that this is due to our choice of t₁ and t₃ values. As given in Fig. 2, t₁ increases the local moment in “case A”, but decreases the local moment in “case B”, while t₃ has completely the opposite effect on the local moments in “case A” and “case B” in comparison with t₁. Our choice of t₁ and t₃ can mimic the magnitude of the parameter in BLG, so we can draw a conclusion that in BLG, “case A” always shows a larger local moment than “case B”. However, in other Bernal-stacked BHLs, the local moment for “case A” and “case B” shall be compared after careful calculations.

We can see from Fig. 3 that the curves for “case A” and “case B” with the same values of parameters intersect. Note that the local moment is determined by the competition between the impurity charge ⟨n_d⟩ and the double occupancy as given in Eq. (5), and the intersection is mainly caused by the crossover of ⟨n_d⟩ curves as μ decreases. In a normal metal, it is well known that the broadening of the impurity energy level depends mainly on the hybridization strength (since the DOS is usually set as a constant). However, in our system where the DOS or the LDOS is not a constant obviously, the broadening of impurity energy level is also affected by the LDOS on the two sublattices, and the renormalized impurity energy level is influenced by the on-site Coulomb repulsion U. As μ is lowered, the impurity occupation number n_d for “case A” and “case B” crosses over each other leading to the intersection point in the local moment, but the intersection has not been observed in the ⟨n_d↑n_d↓⟩ curves.

The results of versus μ for various values of hybridization interaction V are given in Fig. 4. The parameters are chosen as β = 40t⁻¹, U = 1.6t, t₁ = 0.2t, t₃ = 0.1t, and ε_d = −U/2. As in Fig. 4, the lines from top to bottom are listed in the order of increasing V. Circles represent the results for “case A” and triangles for “case B”. We can easily see the same tendency as that given in Fig. 3 when the chemical potential is lowered, and that the local moment decreases monotonically as μ is lowered. The local moment is decreased as V increases. The hybridization interaction mixes the localized states on the impurity site with the itinerant electrons, such that large V lowers the local moment on the magnetic impurity. One need to notice that the local moment for “case A” seems to be larger than that in “case B” when μ is around half-filling for all the choices of V, and this is mainly caused by our choice of t₁ and t₃. In a BHL, if we choose different values of t₁ and t₃, the results for “case A” and “case B” shall change accordingly, and the local moment for these two cases are determined by the competition of the two inter-layer hopping energies.

Fig. 4.

Figure Option ViewDownloadNew Window

Fig. 4. (color online) versus μ for various values of V. The parameters are chosen as β = 40t−1, U = 1.6t, t1 = 0.2t, t3 = 0.1t, and εd −μ = −U/2. The lines from top to bottom are listed in the order of increasing V. Circles represent the results for “case A” and triangles for “case B”. The distinction of for a given value of V is determined by the interlayer hopping energies, and larger local moment is favored for smaller V.

In Fig. 4, we can see similar intersection behavior of the curves for the same values of parameters as those shown in Fig. 3. This is also caused by the intersection of the ⟨n_d⟩ curve as μ decreases. The broadening of the impurity energy level is influenced by the LDOS and hybridization, while the on-site Coulomb repulsion U affects the renormalized impurity energy. As μ is lowered, the impurity occupation number ⟨n_d⟩ for “case A” and “case B” crosses over each other leading to the intersection point in the local moment.

Figure 5 shows the spectral densities A(ω) at 1/T = 40t⁻¹ with various values of U and V. We fix μ = 0, ε_d = −U/2, t₁ = 0.2t, and t₃ = 0.1t, so particle–hole symmetry is preserved and we have A(−ω) = A(ω). It is well known that according to the Hartree–Fock solution of a magnetic impurity in a normal metal, the locations of the peaks are independent of V and the separation of two peaks is about U. However, in Fig. 5 we see that the separations of two side peaks in A(ω) are much smaller than Coulomb repulsion U. This property directly originates from the non-constant DOS in the BHL with Bernal stacking, and is consistent with the results given by the numerical renormalization group study.^[43] When V increases, the A(ω) peaks move toward the charge neutral point. These results are also in agreement with those obtained in single layer graphene.^[36] Near the Dirac point, the spectral densities have the same order as that of the LDOS, because when V is induced, the eigenstates in the host system can greatly hybridize with impurity orbit, and as V increases, the LDOS of the host material influences the impurity orbit more strongly.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) A(ω) with respect to ω for (a) V = 0.75t, U = 1.6t, (b) V = 1.0t, U = 1.6t, and (c) V = 1.0t, U = 2.4t. Here εd = −U/2, β = 1/T = 40t−1, t1 = 0.2t, and t3 = 0.1t in all cases.

In Fig. 5(a) where V = 0.75t and U = 1.6t, we can see that the differences in spectral densities in the three cases are relatively small near ω = 0. In Fig. 5(b), as V is increased from 0.75t to 1.0t, the peaks in A(ω) move towards ω = 0 and thus the A(ω) for “case B” increases dramatically while no such changes are induced in the other two cases. We also study the A(ω) with a larger value of Coulomb repulsion U = 2.4t and V = 1.0t with particle–hole symmetry in Fig. 5(c). It is clear that the separations of two peaks in A(ω) increase but are still much smaller than U. Near the charge neutral point, the spectral densities for the three cases have the same order as that given in Figs. 5(a) and 5(b).

In single layer graphene, the linear dispersion relation and vanishing DOS at the Fermi energy lead to unusual Friedel oscillations^[44] as well as RKKY interactions.^[45,46] Both the Friedel oscillations and RKKY interactions in graphene decay faster than those in normal metals, and our results of charge and spin correlations between the magnetic impurity and the conduction electrons can reflect these behaviors. Especially, the spin–spin correlation shows the spatial distribution of the spin compensation cloud around the magnetic impurity, and has been discussed recently in more complicated systems such as topological semimetals^[47,48] and topological surface states^[49] with spin–orbit couplings.

In Fig. 6, we present the results of charge–charge correlation C_i and spin–spin correlation S_i between the localized state and the conduction electrons, which are given in Eqs. (6) and (7), respectively. We choose the values of the parameters as V = 1.0t, U = 0.8t, ε_d = −U/2, and β = 1/T = 40t⁻¹. μ is fixed at the zero point, so the system preserves particle–hole symmetry. The magnitudes of inter-layer hopping energies are selected as t₁ = 0.2t and t₃ = 0.1t, which are close to the values in a real BLG system. The impurity location is labeled as i = 0, such that the sites with even indices are on the same sublattice as the site where the impurity is added, and those with odd indices are sites on the opposite sublattice.

Fig. 6.

Figure Option ViewDownloadNew Window

Fig. 6. (color online) (a) The charge–charge correlation Ci and (b) the spin–spin correlation Si versus site i. i is along a zigzag direction in layer 1 of the BHL, and the impurity is located on top of site i = 0. Insets (a1) and (a2): details of C0 and C1. Insets (b1) and (b2): details of S0 and S1. In inset (b3), the sites along the red dashed (blue dotted) zigzag lines indicate the sites we consider for “case A” (“case B”), where the blue (red) circles represent the sites in sublattice A (B) in layer 1. The values of the parameters are chosen as V = 1.0t, U = 0.8t, εd = −U/2, and β = 1/T = 40t−1.

Firstly, C_i and S_i for all the three cases are relatively short-ranged and their magnitudes decay rapidly with respect to the location i. This behavior is very similar to those observed in graphene.^[44–46] Secondly, we can see that C_i presented in Fig. 6(a) lacks oscillations. This is because at μ = 0 the BHL is half-filled that the charge exchange between sites is greatly suppressed. In addition, from the spin–spin correlation S_i for μ = 0 given in Fig. 6(b), we can see that the impurity spin is antiferromagnetically correlated with the conduction electron spins on the same sublattice, and ferromagnetically correlated with those on the opposite sublattice. Due to the bipartite nature of the lattice and the localization of impurity spin, the system would show weak spin fluctuation around the magnetic impurity at half-filling.

Since both the charge and spin correlations decay very fast along the zigzag lines, the screening of the conduction electrons on the impurity charge as well as spin is mainly contributed from the sites near the impurity location. If we look at the correlations in details from Fig. 6(a₁), we find that “case B” shows larger on-site correlation |C₀| than “case A”. Similarly, the on-site spin–spin correlation S₀ shown in Fig. 6(b₁) is stronger for “case B”. We can relate the local moment formation for “case A” and “case B” with the short-range correlation strength. The magnitude of the on-site correlation function depends on the LDOS of the carbon sites. Note that the distinction for “case A” and “case B” is determined by our choice of t₁ and t₃. In this case, the dimer sites have lower LDOS than the non-dimer sites, so both C₀ and S₀ have smaller magnitudes in “case A” in comparison with “case B”. Weaker screening leads to larger local moment formed on the impurity site, so at μ = 0 “case A” always exhibits larger local moment than in “case B”, which is observed in both Figs. 3 and 4.

We also focus on the behavior of on-site correlation functions C₀ and S₀ with μ moving below the zero point in Fig. 7. In this case the particle–hole symmetry is broken and the filling is shifted away from one. Here the parameters are chosen as V = 1.0t, U = 0.8t, ε_d = −U/2, and β = 40t⁻¹. As μ is tuned below the Dirac point, the amplitude of on-site charge correlations C₀ in Fig. 7(a) increase because the occupancy of the impurity orbit and conduction band shift from half-filling so their charge exchange enhances. The differences of C₀ for the three cases are the most obvious at μ = 0. As μ is lowered, the differences become smaller and at μ ≈ −0.2t, three curves touch each other. We can see the same behavior of S₀ in Fig. 7(b) as that in C₀. The spin–spin correlations function directly depends on the LDOS^[9] of the conduction band, where t₁ decreases the LDOS on sublattice A and increases the LDOS on sublattice B.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. (color online) (a) The charge–charge correlation C0 and (b) the spin–spin correlation S0 versus μ. Here we use V = 1.0t, U = 0.8t, εd = −U/2, β = 1/T = 40t−1, t1 = 0.2t, and t3 = 0.1t.

We have systematically investigated the asymmetric response of local moment to the inter-layer hoppings in a Bernal-stacked BHL using the non-perturbative quantum Monte Carlo method. Due to the inter-layer hopping energies t₁ and t₃, the two sublattices in the same plane become inequivalent, and this inequivalence is reflected by the asymmetric response of the local moment to the inter-layer hoppings. Interestingly, both t₁ and t₃ show completely opposite effects in the two cases, namely, when the magnetic impurity is coupled to the dimer sites (“case A”) and the non-dimer sites (“case B”), respectively. When the magnetic impurity is coupled to a dimer site, t₁ increases the local moment while t₃ reduces it. Otherwise if the magnetic impurity is located on a non-dimer site, the local moment is reduced as t₁ increases, but is reinforced as t₃ increases. In general, the magnitude of the local moment is determined by the competition between these two inter-layer hopping energies. By varying the model parameters, we can tune the local moment between relatively large and small values, and we find that the asymmetric response persists in a very wide parameter range. The asymmetry between the two sublattices in the same plane in a Bernal-stacked BHL is reflected by the distinctions of LDOS, so it can be measured by the scanning tunneling microscopy (STM) experiments. The spin and charge correlation functions, which also reflect the asymmetric response in short-range, can possibly be measured by the STM and spin-polarized STM.^[50–52] Moreover, the asymmetric response of local moment formation can be explored in experiments through the measurements of magnetic susceptibility. Our results have possible application value to the research of Kondo effects in the class of systems whose electronic structures can be described by a BHL, e.g., the hexagonal boron nitride bilayer (h-BN),^[6] the commensurate h-BN graphene bilayer,^[7] the Mn⁴⁺ ions in Bi₃Mn₄O₁₂(NO₃)^[8] and in cold atom experiments.